\section{Introdução}
 
O sistema sob investigação é o modelo XY clássico ferromagnético apresentado na Seção \ref{sec:xy}  com diluição por sítio do tipo quenched, descrito pelo hamiltoniano
 \begin{eqnarray}
H = -J\sum_{\langle i,j\rangle}\epsilon_i \epsilon_j [S_{i}^{x}S_{j}^{x}+S_{i}^{y}S_{j}^{y}] ,
\label{hamxy}
\end{eqnarray}
%
onde $\vec{S_i}$ representa vetores de spins clássicos tridimensionais  $\vec{S_i}=(S^x_i,S^y_i,S^z_i)$  com  
$S_i^2=(S^x_i)^2+(S^y_i)^2+(S^z_i)^2=1$, a soma é executada sobre todos os pares de primeiros vizinhos em uma rede cúbica simples, e $\epsilon_i$  são variáveis aleatórias,  descorrelacionadas do tipo temperada, que representa a existência de duas classes de partículas no sistema, as partículas magnéticas  ($\epsilon_i=1$) e as não magnéticas ($\epsilon_i=0$)  

As variáveis $\epsilon_i$ obedecem a seguinte distribuição de probabilidade 
%
\begin{equation}
P(\epsilon_i) = p\delta(\epsilon_i-1) + (1-p)\delta(\epsilon_i),
\label{eq:pro}
\end{equation}
%
onde $p$ é a concentração de sítios magnéticos, tal que $p = 1$ corresponde ao sistema puro. 

No Capítulo \ref{cap:cg} já se discutiu,  em linhas gerais, os métodos empregados em simulações Monte Carlo e  análise de dados. Descreve-se com mais detalhes nesta seção o procedimento utilizado neste modelo. 

Para cada configuração de concentração $p$, temperatura $T$ e tamanho de rede são geradas amostras aleatórias $\{\epsilon \}$ usando a distribuição da Eq.(\ref{eq:pro}). Para toda observável $Q$,  foi calculado primeiramente a média térmica de uma amostra $\langle Q_{\{\epsilon\}} \rangle$  e com o resultado das diferentes amostras é calculada a média configuracional com
%
% $\left[ {\left\langle {Q_{\left\{ \epsilon  \right\}} } \right\rangle } \right]_{{\text{av}}}$.
%
%Para encontrar o resultado final, para cada diluição, temperatura e tamanho da rede, o algoritmo estima $\langle Q_{\{\epsilon\}} \rangle$ de quantidades termodinâmicas, para uma dada distribuição aleatória de sítios $\{\epsilon\}$ onde é calculado a media sobre todas as diferentes realizações de desordem 
%                            
\begin{equation}
\left[ {\left\langle {Q_{\left\{ \epsilon  \right\}} } \right\rangle } \right]_{{\text{av}}}  = \frac{1}
{{\# \left\{ \epsilon  \right\}}}\sum\limits_{\left\{ \epsilon  \right\}} {\left\langle {Q_{\left\{ \epsilon  
\right\}} } \right\rangle } ,
\label{av}
\end{equation}
%
onde  $\#\{\epsilon\}$ é o número total de realizações considerado. 

A média sobre as desordens, foi tomada usando amostras independentes. Para todas as amostras são realizadas $1 \times 10^3$ MCS por spin para termalização. E para o calculo da média foi realizado $5 \times 10^4$ MCS. 
O tamanho da rede varia de $L=10$, $20$, $30$, $40$, $50$, estes valores foram escolhidos para que $p \times L^3$ fosse um número inteiro. 


Para cada amostra de uma dada configuração de diluição, foi utilizado um algoritmo de Monte Carlo Híbrido consistindo de três passos do algoritmo metropolis (M) cinco passos de uma versão não ergódiga do algoritmo de Wolff (W)~\cite{Wolff1989,newman1999monte}, e dois passos do algoritmo de superrelaxação a energia configuracional constante (O)~\cite{creutz1987overrelaxation, pawig1998monte}.  Os passos individuais são misturados automaticamente no programa e então uma seqüência é gerada, como por exemplo (MWWMOWMWWO), este é considerado então um passo de Monte Carlo Híbrido. O algoritmo de  Monte Carlo Híbrido foi implementado e mostrou redução nas correlações entre as sucessivas configurações de spin na simulação. Próximo da temperatura de transição também foi utilizado a técnica do histograma simples para obter as correspondentes quantidades Termodinâmicas.  
Inicialmente foi calculado a susceptibilidade no plano e o cumulante de Binder 
%
\begin{eqnarray}
m_{xy}&=&\frac{1}{L^{2}}\sum_{i}^{L^{3}}[(S_{i}^{x})^{2}+(S_{i}^{y})^{2}],\\
\label{mag}
%
\chi_{} &=& L^{2} \frac{\langle m_{xy}^{2}\rangle -\langle m_{xy}\rangle ^{2}}{T},\\
\label{chi}
%
u_{4}&=&1-\frac{\langle m_{x}^{4}\rangle}{3\langle m_{x}^{2}\rangle ^2},
\label{U1}
\end{eqnarray}
%
onde $L$ é o tamanho linear do sistema da rede cúbica estudada, $T$ é a temperatura dada em unidades de $J/k_B$, sendo $k_B$ a constante de  Boltzmann, e para o cumulante tem-se $m_\alpha^m = \left({\frac{1}{L^3}}\sum_{i=1}^{L^3}{S_i^\alpha}\right)^m$, onde $\alpha=x,y,$  ou $z$ (a componente $x$ tem se mostrado ser mais apropriada para $u_4$ \cite{Landau1999}).

O expoente crítico foi estimado utilizando a teoria de escala de tamanho finito sobre a média do modulo do parâmetro de ordem : $\langle \sqrt{m_{x}^{2} + m_{y}^{2}} \rangle$ , sobre a média do quadrado do parâmetro de ordem $\langle m_{x}^{2} + m_{y}^{2} \rangle$, e na derivada da temperatura deste último. Em $T=T_c$ é encontrado as relações já discutidas na Seção \ref{sec:escala}, 
%
\begin{eqnarray}
 X_1&=& L^{ - 3} \left\langle {m_x^2  + m_y^2 } \right\rangle  = X_{10} L^{  \gamma /\nu } 
 (1+X_{11}L^{-\omega}) 
\label{gamma} \\
%
X_2&=&  L^{ - 3} \left\langle {\sqrt {m_x^2  + m_y^2 } } \right\rangle  = X_{20} L^{ - \beta /\nu }
(1+X_{21}L^{-\omega})
\label{beta}\\
%
X_3&=& \frac{\partial }{{\partial T}}\ln \left\langle {m_x^2  + m_y^2 } \right\rangle = X_{30} 
L^{  1/\nu }(1+X_{31}L^{-\omega}),
\label{nu}
\end{eqnarray}
%
onde $\beta$, $\gamma$, e $\nu$ são os expoentes críticos, $X_{ij}$ é uma constante não universal e $\omega$ é o expoente de correção de escala.



\section{Resultados}
% \subsection{Rede cúbica}  
%\subsection{Quantidades termodinâmicas}
Sabe-se que simulação computacional de sistemas desordenado, próximo da temperatura de transição , apresenta pequenas variações nos valores das quantidades termodinâmicas para diferentes amostras realizadas. Como exemplo, na Figura \ref{xxsL20p85} está apresentado os valores da susceptibilidade para cada uma das amostras em uma temperatura próxima da temperatura crítica, para concentrações de $p=0.95$ e $ p=0.85$,
e tamanho de rede de $L=20$. Na figura os pontos são os valores da susceptibilidade para a $n$-ésima amostra, a linha cheia é a média de $n$ amostras. Pode-se observar que a distribuição dos pontos em torno da média é quase simétrica e torna-se mais disperso a medida que aumenta a concentração de átomos não magnéticos. Na Figura \ref{disL20p85} observa-se que a distribuição dos valores da susceptibilidade para diferentes amostras se aproxima de uma Gaussiana e que a média e a mediana estão muito próximas. Assumindo a aproximação gaussiana para a distribuição, pode-se calcular o intervalo de confiança \cite{cardy1996scaling}.  Com esta análise percebe-se que a média é estável acima de $100$ realizações de desordem e que o intervalo de confiança diminui lentamente com o número de amostra. Com $ 100 $ amostras o intervalo de confiança fica abaixo $ 0.6 \%$ para $ p =  0.85$. Uma vez que o custo computacional para reduzir esse erro é muito alto, especialmente para redes maiores, foi utilizado neste trabalho $100$ amostras para todas as configurações.


%
\begin{figure}[htb]
\begin{center}
\includegraphics[width = 10cm]{figuras/xydiluido/fig1.jpg}
\end{center}
\caption{Susceptibilidade magnética $\chi_{\{\epsilon\}}$ em diferentes amostras,  para  
$p=0.95$ e $p = 0.85$ e  $L = 20$, obtida na temperatura próxima do ponto onde a média  $[\chi_{\{\epsilon\}}(T)]_{av}$  tem um máximo. A linha sólida é a média de acordo com a equação. (\ref{av}).}
\label{xxsL20p85}
\end{figure}
%
%
\begin{figure}[htb]
\begin{center}
\includegraphics[width = 8cm]{figuras/xydiluido/p085.jpg}
\end{center}
\caption{Distribuição da susceptibilidade magnética $\chi_{\{\epsilon\}}$  para  
$p=0.95$ e $p = 0.85$ e  $L = 20$.  A linha sólida é um ajuste por uma gaussiana.}
\label{disL20p85}
\end{figure}
%
%\subsection{Expoentes críticos}
  
A fim de estudar o comportamento do expoente crítico com diluição e boa precisão, restringiu-se o estudo a concentrações maior que $0.85$, pelo fato do intervalo de confiança decrescer com o aumento da concentração de sítios magnéticos, como pode ser vistos na Figura \ref{xxsL20p85}. Apesar da restrição no intervalo de concentração é possível, a partir desse, obter conclusões gerais para o modelo.   
Como exemplo típico, nas Figuras \ref{cx1-3} apresenta-se o resultado obtido para $p=0.97$ para o cumulante, e as quantidades termodinâmicas dadas pelas Equações (\ref{gamma})-(\ref{nu}).  Tomando o cruzamento das duas maiores redes no cumulante, encontra-se $T_c^{u_4}=1.5008(1)$. 

Pode-se, no entanto, a partir da escala de temperatura da Figura \ref{cx1-3} perceber que ainda há uma dependência do tamanho da rede no $u_4$. Contudo usando uma estimativa de para o expoente $\nu$ e ajustando os dados para uma dependência de escala de tamanho finito (sem correção de escala) dos cruzamentos das redes $L=30,40,50$  com $L=20$ foi encontrado  $T_c^{u_4}=1.5009(2)$. Isso mostra que a estimativa acaba por não ser tão diferente da anterior. Logo, para obter a temperatura crítica sem a necessidade de qualquer expoente usou-se apenas a temperatura crítica do cruzamento de $u_4$.
  
Na figura \ref{cx1-3}b-d mostra um exemplo do comportamento das quantidades $X_1$, $X_2$, e $X_3$ com o tamanho da rede na temperatura $T_c^{u_4}=1.5008$. Dessas, foram estimados os expoentes críticos ,inicialmente negligenciando o termo de correção de escala de tamanho finito, i.e. $X_{11}=X_{21}=X_{31}=0$ nas Equações. (\ref{gamma})-(\ref{nu}). Os resultados são mostrados na Figura \ref{cx1-3} e os expoentes estão apresentados na tabela \ref{tab}. Pode-se observar que os pontos se ajustam bem a uma linha reta, o que implica que correções de escala não são importantes para o tamanho de rede utilizado. Usando o mesmo procedimento foram obtidos os expoentes críticos para outras concentrações. Os respectivos resultados estão apresentados também  na Tabela \ref{tab}. Na Tabela também estão inclusos os resultados obtidos para $p=1$ que foram comparados com os da referencia \cite{PhysRevB.74.144506}. Pode-se observar que há um completo acordo com os expoentes do modelo puro. Isto é consistente como o critério de Harris, apesar do expoente do calor específico ser bastante pequeno, mas é de fato negativo. O Cumulante universal de  Binder também concorda com este critério. 
 %
\begin{figure}[htb]
\begin{center}
\includegraphics[width = 10cm]{figuras/xydiluido/cx1-3.jpg}
\end{center}
\caption{Quantidades termodinâmicas para o modelo XY diluído como $p=0.97$. (a) Cumulante de Binder como função da temperatura para diferentes tamanhos de rede. (b)-(d) Escala de tamanho finito de $X_1$, $X_2$, 
e $X_3$, respectivamente, com o correspondente ajuste  (linha cheia)  sem as correções de  escala.  As barras de erro são menores que o tamanho do símbolo. }
\label{cx1-3}
\end{figure}



\begin{table}[htbp]
\caption{Expoentes críticos do modelo XY com diluição por sítios para diferentes concentrações. Também é dado cumulante de Binder e a temperatura crítica obtida do máximo da susceptibilidade e do cumulante de Binder. 
%The exponents for the pure model are from reference \cite{massimo} and
%the cumulant from reference \cite{prb60}.
}
\begin{center}
\begin{tabular}{llllll}
\hline \hline
 & $p=1$\cite{PhysRevB.74.144506}& $p=1$        & $p=0.97$ & $p=0.95$ \\ \hline\hline
 
 $\alpha$ & -0.0151(3) & -0.0037(104)  & -0.0134(38) & -0.0055(106)\\
% & -0.0151(3) & -0.041(30)  & -0.009(63) & -0.063(48) \\ \hline

$\beta$~~~ & 0.3486(1) &0.3423(22) & 0.3443(13)~~ & 0.3460(36)~~~~ \\
% ~~~ &  0.3486(1)~~ &0.344(25) & 0.351(27)~~ & 0.356(32)~~~~ \\ \hline

$\gamma$~~~ & 1.3178(2) & 1.3164(75)  & 1.3194(35) & 1.314(10)\\
%~~~ & 1.3178(2) & 1.338(15)  & 1.318(15) & 1.323(16)\\ \hline

$\nu$ & 0.6717(1) & 0.6679(35)      & 0.6711(13) & 0.6685(35)\\
% & 0.6717(1) & 0.68(1)      & 0.670(21) & 0.687(18) \\ \hline

%$\omega$ &  & 0.790(11)  & 0.807(29) & 0.80(2) \\ \hline

$T_c^{u_4}$&~~~& 1.55177(9)&1.5008(1)&1.46631(31) \\ \hline
$T_c^{X_1}$&~~~ & 1.55184(7)& 1.50060(23)& 1.46615(29)\\ \hline
$u_4$ & 0.3789(15)\cite{PhysRevB.60.3375}&0.3808(10) & 0.3764(15) &0.3798(64) \\ \hline\hline

\end{tabular}
\end{center}
\label{tab}
\end{table}

Uma vez que foram calculados os expoentes críticos pode ser obtida uma estimativa adicional mais precisa da temperatura crítica do máximo da susceptibilidade. A localização dos picos define uma temperatura crítica (efetiva) dependente do tamanho da rede $T_c(L)$  que escala de acordo a equação
%
\begin{equation}
T_c (L) = T_c  + \lambda L^{-1/\nu},
\label{eq:tc}
\end{equation}
%
onde o expoente crítico universal $\nu$ foi obtido das Equações (\ref{gamma})-(\ref{nu}) e a constante não universal Fig. \ref{fittc}  mostra o ajuste dos pontos de acordo com a relação de escala dada acima Eq.(\ref{eq:tc}). O correspondente valor da temperatura crítica extrapolado é dado também na Tabela \ref{tab}. Nota-se que a temperatura crítica do máximo da susceptibilidade e do cumulante de Binder são bastante semelhantes. Assim, em alguns casos a temperatura crítica só foi determinada pelo cruzamento do cumulante.   Também se poderia tentar analisar as correções de escalas 
na equação. (\ref{gamma}) - (\ref{nu}) para este caso. Contudo os termos de correção são muito pequenos, o que tende a dificultar a determinação de $ \omega $ com uma boa precisão.  Em todo caso foi obtida uma estimativa de $\omega\sim 0.80(10)$, que é comparável ao $\omega = 0.785(20) $ da referência \cite{PhysRevB.74.144506}.

\begin{figure}[htb]
\begin{center}
\includegraphics[width = 10cm]{figuras/xydiluido/tc.jpg}
\end{center}
\caption{ Temperatura $T_c(L)$ como função de tamanho linear da rede $L$. A linha sólida mostra o ajuste dos dados de acordo com a relação de escala dada na Equação (\ref{eq:tc}). Os mesmos valores estão apresentados na Tabela \ref{tab}. As barras de erros são menores que o tamanho dos símbolos.}
\label{fittc}
\end{figure}
 
O diagrama de fase obtido está apresentado na Figura \ref{dia}. Pode-se observar que o erro em $T_c$ cresce com a diluição. Este fato está relacionado à necessidade de um alto número de amostras para sistemas altamente diluídos pois muitas configurações não fazem um agrupamento de spins percolado.  
A redução inicial da temperatura $T_c(p)$ com a diluição pode ser definida por $S=\left. {\frac{d}{{dp}}\left( {\frac{{T_c (p)}}{{T_c (1)}}} \right)} \right|_{p = 1}$. 


O valor de $S$ encontrado neste trabalho ($ S=1.0965(39)$) é semelhante ao $S=1.202$ calculado por aproximação de campo efetivo \cite{tucker1996effective} na versão quântica de spin-1/2 do modelo. Estes valores são bastante diferentes do resultado experimental onde foi encontrado $ S=0.24(2)$ \cite{defotis2008dependence}. Contudo para alta diluição não magnética ha uma boa concordância na inclinação correspondente. O presente comportamento experimental sugere que um acoplamento extra de interação de supertroca possa ser a causa da baixa inclinação inicial com a diluição, como já foi discutido na referência \cite{defotis2008dependence}. Isto é similar ao que acontece nas ligas Fe-Al, onde  interação de supertroca explica a inclinação menos brusca em baixas concentrações de Al \cite{PhysRevB.61.3188,dias2009ising} 

\begin{figure}[htb]
\begin{center}
\includegraphics[width = 10cm]{figuras/xydiluido/diag+exp.jpg}
\end{center}
\caption{Diagrama de fase do modelo 3D XY com diluição por sítios do tipo quenched no plano temperatura reduzida $T_c(p)/T_c(1)$ versus concentração de sítios magnéticos   $p$. O quadrado marca o local do ponto de percolação $p_c$\cite{stauffer1994introduction}. A linha pontilhada é uma guia para os olhos. O destaque mostra a inclinação  para baixa concentração, comparado com resultados experimentais da referência\cite{defotis2008dependence}. A linha pontilhada no caso, dá a correspondente inclinação na região indicada. }
\label{dia}
\end{figure}

%\subsection{Rede quadrada}



